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Abstract 

We analyze the finite volume spatial convergence combined with Exponential 
Time Differencing of order one (ETD1) for temporal discretization of a non 
linear Advection-Diffusion-reaction equation (ADR) modeling transport in 
porous media. We derive the L 2 estimate under the assumption that the non 
linear reaction is globally Lipschitz. We illustrate the theoretical proof by 
some applications in 2D and 3D highly anisotropic and heterogeneous porous 
media including the SPE10 case pQ. We compute the exponential of the non 
diagonal matrices arriving in the finite volume spatial discretization with the 
Fast Leja points technique and the Krylov subspace technique. 
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1. Introduction 

Flow and transport are fundamental in many geo-engineering applica- 
tions, including oil and gas recovery from hydrocarbon reservoirs, ground- 
water contamination and sustainable use of groundwater resources, storing 
greenhouse gases (e.g. CO2) or radioactive waste in the subsurface, or mining 
heat from geothermal reservoirs. In porous media, a single phase transport 
process is described by the equation 
dX 

^(x)- = V-(D(x)VI)-V-(q(x)I) + fl(t,x,I) (x,t) Gfix[0,T],(l) 
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where Q is an open domain of IR d , d G {1,2,3}, is the porosity (void 
fraction) of the rock, and D is the symmetric dispersion tensor, X is the 
unknown concentration of the contaminant, q the Darcy's velocity and R 
the reaction and source term. 

Finite element, finite volume or the combination finite element-finite volume 
methods are mostly used for space discretization of the problem ([!]) while ex- 
plicit Euler, semi implicit and fully implicit methods are usually used for time 
discretization [51 [7J [91 [10]. Due to time step size constraints, fully implicit 
scheme is more popular for time discretization for quite a long time com- 
pared to explicit Euler and semi implicit methods. This method, however, 
need at each time step a solution of large systems of nonlinear equations. 
This can be the bottleneck in computations. In recent years, exponential in- 
tegrators have become an attractive alternative in many evolutions equations 
[21 El [H [191 El HU H21 HZl I2QI HH H2l [H]. In contrast to classical methods, 
they do not require the solution of large linear systems. Instead they make 
explicit use of the matrix exponential and related matrix functions. The 
price to pay is the computing of the matrix exponential functions of the 
non diagonal matrices, which has revived interest and significance progresses 
nowadays [H HS1 HB H2J [201 HB [Lj . 

In this work, we combine a finite volume method with the first order 
exponential time differencing scheme of order 1 (ETD1). Although both 
discretization techniques have been together used for solving evolutionary 
problems like Q (see [21 El i]), a proper combination of rigorous proof of 
them has been lacking so far. 

The paper is organised as follows. In Section [2] we present the semi 
group formulation. In Section [3] and Section [4] we present the space and time 
discretization. We then state and discuss our main result in Section [5] and 
present the proof of the convergence theorem. We end by presenting some 
simulations in Section [6j these are applied both to a linear example where 
we can compute an exact solution as well as a more realistic model (SPE10 
case [T]). 

2. Semi group formulation and well posedness 

Let us start by presenting briefly the notation of the main function spaces 
and norms that we use in this paper. We denote by || • || the norm associated to 
the inner product (•, •) of the Hilbert space H = L 2 (Q). For a Banach space V 
we denote by || • || v the norm of V and L(V) the set of bounded linear mapping 
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from V to V to C. For the sake of simplicity, without loss of generality 
we use the homogeneous Dirichlet boundary condition and assume that the 
porosity <p is constant. During our simulation, we will also investigate the 
general boundary conditions and non constant porosity (p. We assume that 
the Darcy velocity q is known, and satisfies the following mass conservation 
for incompressible fluids V ■ q = 0. The model problem ([!]) is reformulated 
as: find the concentration X(t) G H 1 ^) such that 

{dX/dt + AX = i?(x, t, X) (x, t) G Q x [0, T] 

X(x,0) = X xGfi (2) 

x(x,t) = o (x, t) g on x [0, T] 

where 

AX = A(pc, X) = - V ■ (DVX) + V ■ (q(x)X) 



( J=l J 1=1 



For well posedness of Q, we assume that D is symmetric, D iy j G L°°(f2), gj G 
Z/°°(fi) and there exists a positive constant C\ > such that 

^Aj(x)^j > Ci|e| 2 V£ G lR d xgH Cl >0. (3) 

Set V = Hq(Q), the bilinear form associated to the operator A is given by 

= JA it, D * |r ^ + E^ ) dx u > ve v - ^ 



According to Garding' s inequality (see [191 EI]); there exists two positive 
constants cq and Ao such that 

a(v,v) + c 1 1 ^ 1 1 2 > ^o||^||lfi(r2) Vv e V. (5) 
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By adding cqX in both side of the first equation of ([T| we have a new operator 
that we still call A corresponding to the new bilinear form that we still call 
a such that the following coercivity property holds 

a(v,v)> Ao|M|tfi (n) VveV. (6) 

We will still call the right hand side of the first equation of (|2j) R. 
By Green's formula we have 

a(u, v) = (Au, v) Vu G H] n H 2 {tt) = V(A), \/v G V, (7) 

where T>(A) the domain of the operator A. Therefore the weak form of ^ 
is to find the function X(t) G T>(A) such that 

'(X t ,x) + (AX,x) = (R(X),x) V X GV, tE[0,T] 

X(t) =x . u 



The V— ellipticity (J6J) implies that —A is a sectorial on L 2 {Vt) (see [23| 121]) 
i.e. there exists Ci, 9 G (|tt, 7r) such that 



IKAJ + A)- 1 !!^))^^ AG5 e , (9) 

where 5 e = {A G C : A = pe^, p > 0, < \<j>\ < 9}. 

Then —A is the infinitesimal generator of bounded analytic semigroups 
S(t) := e~ tA on L 2 (fi) such that 

S'(t) :=e~ iA = — [ e tx {\I + A)- 1 d\ t>0 (10) 
2m J c 

where C denotes a path that surrounds the spectrum of —A. 

The coercivity property in ^ implies also that the set of the real part of 
the spectrum of A is non negative, which allows the definition of the fractional 
power of A as: for any a > 

[a) (11) 

where r(a) is the Gamma function of a [23]. We denote by ||.|| Q := ||y4 Q / 2 .|| 
the norm of the space T>(A a l 2 ). 

For the nonlinear reaction term we make the following assumption 
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Assumption 2.1. [Lipschitz condition for R] 

The nonlinear R is continuous and differentiable respect to X with respect 
to the variable X, with bounded differential. Therefore there exists a positive 
constant L > such that 

\R(x,t,u) - R(x,t,v)\ < L\u-v\ \/u,veR, x £ H, t e [0,T] 

(12) 

|| R(Y ) -R(Z) 1 1 < L\\Y- Z\\ VY,Z e L 2 (Q), 

where with a slight abuse of notation, R denotes the nonlinear operator 
X^R(;;X). 

By Duhamel's principle we may represent solutions of ^ by the following 
integral equation 

X(t) = S(t)X + I S{t- s)R{s,X{s))ds, te[0,T]. (13) 
Jo 

Applying the contraction mapping principle in the topology of the Banach 
space C([0, T], H 1 ^)) to the integral equation 
[221 Theorem 6.3.1] ensure the existence and uniqueness of X. 

3. Finite volume for space discretization 

A cell-centred finite volume methods for heterogeneous and anisotropic 
diffusion problems remains a challenging problem. An active area of research 
is to make the approximation of the diffusion flux more efficient and simple 
as possible (see [21] for the references). The finite volume method is widely 
applied when the differential equations are in divergence form. To obtain 
a finite volume discretization, the domain Q is subdivided into subdomains 
(A)ieZ) % being the corresponding set of indices, called control volumes 
or control domains such that the collection of all those subdomains forms 
a partition of Q. The common feature of all finite volume methods is to 
integrate the equation over each control volume A4, i G X and apply Gauss's 
divergence theorem to convert the volume integral to a surface integral. For 
our parabolic problem (|2]), finite volume methods differ in the way they 
approximate the diffusion flux F = — DVX. Two techniques are mostly used: 
the finite volume with two-point flux approximation (TPEA) (see [TJ [21]) and 
the finite volume with multi-point flux approximations (MPFA)( [25l 126] ). 



13) [23, Theorem 3.3.3] or 
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An advantage of the two-point approximation is that it provides mono- 
tonicity properties, under the form of a local maximum principle. It is ef- 
ficient and mostly used in industrial simulations. In this paper we use the 
TPFA as developed in [7J section 11 page 815]. The main drawback of TPFA 
is that it is applicable in the so called " admissible mesh " and not in general 
mesh. 

Definition 3.1. [An admissible mesh [1]} 

An admissible mesh T for problem with the full diffusion tensor D is 
given by: 

• A set {Ai\ i( - X of control volumes such that Q = U with the corre- 
sponding local inner product induced by D^ 1 where 

D ^msk)Z, DW<iX - 

• The corresponding set of center points {xj} ie2 - such that 

(a) Xj G Ai, i 6 I. 

(b) Xj is the intersection of the straight lines perpendicular to the 
boundary of A^ with respect to the inner product induced by D^ 1 . 

Let h be the maximum mesh size of T ■ We denote by l~h a dual Delaunay 
triangulation of T i.e. the Delaunay triangulation where {xj} ig2 - is the set of 
vertices. 

Let us illustre Definition 13.11 to make it more understandable. 

Example 1. • In the case where the diffusion tensor D is diagonal and Q 
is a rectangular or parallelepiped domain, any rectangular grid (d = 2) 
or parallelepiped grid (d = 3) is an admissible mesh. The set {xj} is the 
set of centers of gravity of the rectangular grid or parallelepiped grid. 
The inner product induced locally by D^ 1 is equivalent to the standard 
inner product corresponding to the Euclidean norm \.\. This mesh will 
yield a 5-point scheme (d = 2) and 7-point scheme (d = 3) for our 
model problem p|). 

• If d = 2, for isotropic and heterogeneous media (D(x) = 6(x)/ 2 x£fl, 
I2 being the identity matrix of dimension 2) we can define a triangular 
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admissible mesh T to be a family of open triangular disjoint subsets of 
fl such that two triangles having a common edge have also two com- 
mon vertices. The angles of the triangles are assumed to be less than | 
to allow the orthogonal bisectors to intersect inside each triangle, thus 
naturally defining the center point Xj of the control volume A^. The 
finite volume scheme defined on such mesh will yield a J^-point scheme 
for our model problem |I|). The inner product induced locally by D^ 1 
is equivalence to the standard inner product corresponding to the Eu- 
clidean norm |.|. 

To make notation easier, we will identify T to X, then to say Ai G T we 
will say i G T. 

Consider the modified model problem of ^ where cqX is added on both 
sides of the first equation of problem ([2]), Cq is defined in (J5|. Consider an 



admissible mesh T in the sense of Definition |3.1[ Denote by £ the set of 
edges of control volume of T, S int the set of interior edges of control volume 
of T, Xi{t) the approximation of X at time t at the center (or at any point) 
of the control volume j e T and X a {t) the approximation of X at time t at 
the center (or at any point) of the edge a G £. For a control volume i G T, 
denote by Si the set of edges of i, mes(i) the Lebesgue measure of the control 
volume i G T. 

As in |27], integration over any control volume i £ T, using the di- 
vergence theorem to convert the integral over i to a surface integral, finite 
differences for the diffusion flux approximation |7J and the upwind technique 
for the advection flux approximation yields 



mes(z)^-p + ( F iAt) + Qi,aX a>+ (t)) + c mes{t)X t {t) 



D ia = \T>im !a \, Di = — -T- r.D(x)dx, 

mesu 



F iitT [t) = mes(a) D itU - 



mes(i) R(3q,,t,Xi(t)), 

(14) 



%o = " n i,adcr \/i G T, Ver G £j. 

Here nj )(T is the normal unit vector to a outward to i, mes(<r) is the Lebesgue 
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measure of the edge a G Si and d ijCT the distance between the center of i and 
the edge a. 

Since the flux is continuous at the interface of two control volumes % and 
j (denoted by i \ j) we therefore have F ijCT (t) = —Fj a {t) for a = i \ j, which 
yields 



FUt) = -ra (Xj(t) - Xi(t)) = ( Xj (t) - Xi(t)) , a = i | j 



d 



r a = mes(cx) 



Di y(T F)j y(T 



h3 



Di,a di^ a -\- Dj a dj :Cr 



(transmissibility through a) 



(15) 



with 



\i a = d 



Di,aDj,cr 



Di t a di^ a -\- Dj a dj t<T 



(16) 



where dij is the distance between the center of i and center of j. We will also 
denote by d a the distance d^ or d ija for a — i \ j or a — Si fl dfl respectively. 
Notice that for a C d£l, we can also write 

F ita (t) = -TriXjW-Xitt)) 

mes(a)/v 



di,a 



(Xjit) - Xi(t)) 



with 



{ X j (t)=X a (t)=0 
mes(cr)A, CT 



Trr = 



d. 



1,0 



The upwind term for advection flux X a>+ is defined as 

Xi{t) if q^^O 



(17) 



Xj{t) if q i , a <0 
Xi(t) if q^^O 

X a (t) if q ha <0 



for a — i I j 



for cr e Si n 0ft. 



(18) 
(19) 



We can write X ai + as 

= r a X t (t) + (1 - rjXfi), a = i\j (20) 

where r a = -(sign(g i|(T ) + 1). 

The finite volume space discretization for the model problem (pj is given 

by 

' ,.,dX i (t) ( mes(a)/i CT 
mes(z)— — - + E ~ A {Xj(t) - Xi(t)) 

al creSi V a i,j 

(21) 

+q ha (r a Xi(t) + (1 - rjXjit))) + c mes(i)Xi(t) = mes(z)#(x 4 , t, X 4 (t)) v ; 
k Xj (t) = 0, = if a c Vz G T. 



The scheme (21 ) clearly indicates the affinity of the finite volume method 
to the finite difference method. However, for the subsequent analysis it is 
more convenient to rewrite scheme (21) in a discrete variational form. 

Multiplying the first equation of (21) by arbitrary numbers Vi G R and 
summing the results over all control volume in T yields 

E mes(i)— — + £ ^ (^iW - x j(t)) 

ieT l aZ creSi V a i,j 

< + %a (r a X t (t) + (1 - r a )X,(t)))] Vi = £mes(z) R(Xi(t))vi, ( 22 ) 

ieT 

k R(Xi(t)) :=R(^,t,X t (t)). 

Let Vh G V denote the space of continuous functions that are piecewise linear 
over the Delaunay triangulation Th (dual of T), then the values Xi(t) and 
can be interpolated in Vh- There are unique functions Xh(t),Vh G Vh such 
that Xh(t)(x.i) = Xi(t) and Uft(xi) = for all ieT, where x* is a center of 
the control volume i G T (x* is also a vertex in Th)- 
Denote by the bilinear form defined by 

( Hies(cr) yUo- .\ 

a h {u h , v h ) = 2^ 1^ [ \ u j ~ u i) + Qi,<r ( r oUi + (1 - r a )u-j) 1 ^ 

ieT a-eSi V / ^23) 

+c mes(i)uiVi \/u h ,v h eV h , 



and by (., .) 0i h the scalar product on C(fi) D Vh denned by 
(u,v) ,h = y^mes(i)ujVi, Ui = u(xj), v { = f(xj), u,v G C(Q). (24) 



- 2 , 



The corresponding norm is the discrete L (Q) norm denoted by ||.||o,/i- It is 
proved in [27] that \\-\\o,h is equivalent to the L (Q) norm |.| in Vh when the 
mesh T is regular. 

Previous results allow us to write the following variational form of our finite 



volume scheme (22) 



(^X h ,<p) 0th + a h (X h (t),<p) = (R{X h (t)),<p) fi, V<peV h , £g(0,T], 
dt (25) 

Xh{$) = Xho, 
with 

R(X h )(t)(^) = R(X t (t)) := R(xi,t,Xi(t)) = R(^,t,X h (t)(^)), Vi G T. 
Consider the operator Ah : Vh Vh such that 

(A h ip, x)o,h = ah(ip, X) W>, X e V^h. (26) 
The semi discrete solution in 14 is then given by: find Xh(t) G VJ, such that 

+ A h X h = P h R(Xh) te(0,T\ 
dt ^ 

X h (0) = x oh 

where Ph is the orthogonal projection defined from u G C(O) to Vh by 

(PhU,x)o,fc=<u,X>o,h VxGVfe. (28) 

4. Exponential integrator for time discretization 



In order to give the corresponding mild form of (27) and build the expo- 



nential integrator scheme let us define the discrete Hq(Q) norm. 
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Definition 4.1. [Discrete Hq(CI) norm /W/ 

Let T be an admissible finite volume mesh in the sense of Definition 3J 



Let X{T) be the space of the functions constant in each control volume ofT. 
For u G X(T), the discrete Hq(Q) norm of u is defined by 

1/2 



\ u \\l,T 



(29) 



where 



D a u 
D„u 



mes(a) 



I u i — Uj | 
I Mi I 



if a = i\j e £ int 
if a G <90. 



During our analysis we make the following assumption as in section 
11 page 815]. 

Assumption 4.2. [Regularity of D, q and the admissible mesh T] 

We assume that~D is bounded, the restriction of D to anyi G T belongs to C 1 (z, 
C 1 (Q) and that there exists Ci > and ( 2 > such that 

' Cih 2 < mes(i) < ( 2 h 2 , V* G T, 



(ih < mes(cr) < ( 2 h, Vo" G £ 
k Cih <d a < ( 2 h, Va G £. 



(30) 



Assumption 4.2 allows the following Vh— ellipticity of a^. 

Theorem 4.3. Under the regularity of the admissible mesh T in Assump- 
tion 4j2, there exists a constant a > such that 

ah(vh,v h ) > a |K|li,r e v h- 

Proof. Let b\, b\ and b\ the bilinear forms defined in 14 x V/, by 

mes(cx) \x a 



(31) 



bl(u h ,v h ) 
bl(u h ,v h ) 
b 3 h (u hj v h ) 



EE 



d%, 



(Uj - Ui) Vi 



(32) 



c o y^mes(z) UiVi. (34) 
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Using Assumption 4.2 mainly the regularity of T and the fact that the 
coefficients of the diffusion tensor D are bounded, there exists two constants 
C B (n, Ci, Ca, D) and C' 5 (Q, &, (2, D) such that 



C 5 < /ia = d 



I -J 



^i,a di p -\- Di udi fj 



<C' 5 , a = i\j, 



and 



so that 



c 5 <fx a = Di, ff < c' s , a e £ t n an, 

G> < n„ < C 5 , We £, 



(35) 



(36) 



(37) 



where \x a is defined in (17) and (16). 



Using the fact that the transmissibility given in (15) is symmetric, i.e. 



Ti\j = Tju and reorganizing the summation, we therefore have 



CAK\\lr<^v h )<C' h \\vJ x 



r 



(38) 



Let use some important results from [7]. Indeed as in [7] reordering the 
summation over the set of edges yields 



b 2 h(v h ,v h ) = ^q a (v aj+ - v a>+ 



(39) 



T<=£ 



where 



q* 



Vi if qi, a < 

Vj (or v a ) if q La > 
q ■ n i>a da\. 



a e£ int (oi a e£ l ndn), (40) 

(41) 



Since 



we therefore have 



« + - <_) = W / q • n ^ rf(T ) v *= I"?- q(x) 4 (x)cfc = 0, (43) 
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since we have assumed divergence-free flow. Then 

b 2 h (v h ,v h )>0. (44) 

We also have 

b 3 h (v h ,v h ) = Co\\v h \\ 0ih > 0, (45) 

therefore 

a h (v h ,v h ) > C 5 \\v h \\\,T V v h e V h . (46) 

□ 

The following Vh— ellipticity of a% implies that — Ah is a sectorial on 
L 2 (Q) (uniformly in h) i.e. there exists C\, G (|tt, 7r), such that 

||(AJ + A ft )- 1 || i(L2(fi)) < AG 5 fl , (47) 

where S e = {A G C : A = pe** p > 0, < \<j>\ < 9}. 

The discrete operator — A^ therefore is the infinitesimal generator of 
bounded analytic semigroup (or exponential operator) Sh(t) := e~ tAh on 
Vh such that 

S h (t) := e~ tAh = [ e tx (XI + A^dX , t > (48) 

where C denotes a path that surrounds the spectrum of — A h . As for the 
continuous case, Duhamel's principle implies that the solution of (27) is rep- 
resented by the following integral equations (mild form) 

X h (t) = S h (t)X oh + f S h (t - s)P h R(X h (s))ds, t G [0, T] . (49) 

J o 

For simplicity we consider a constant time-step At > 0. At time t m = mAt G 
[0, T], the mild solution (49) is given by 

X h (t m ) = S h (t m )X oh + / S h (t m - s)P h R(X h (s))ds. (50) 

Jo 
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Then, given the solution at the time t m , we can construct the correspond- 
ing solution at t m +i as 

i-At 

X h (t m+1 ) = S h (At)X h (t rn ) + / S h (At - s)P h R(X h (t m + s))ds. (51) 

Jo 



Note that the expression in (51 ) is still an exact form of X^. The idea behind 
exponential time differencing is to approximate PhR(Xh(t m + s)) by a suitable 
polynomial [13], [H]. We consider the simplest case where PhR(Xh(t m + s)) is 
approximated by the constant PhR(Xh(t m )) and the corresponding scheme 
(ETD1) is given by 

X n+1 = e -AtA hX n + At ^ AtAh )P h R(X^,t m ) (52) 

where 

1 r At 

ipx(-AtA h ) = (-At A,)- 1 (e~ A ^ - J) = — J e^ At ^ds. 



Note that the ETD1 scheme in (52) can be rewritten as 

X^ +l = X™ + At^(-AtA h )(-A h XJ? + P h R(XJ?)). (53) 

This new expression has the advantage that it is computationally more effi- 
cient as only one matrix exponential function needs to be evaluated at each 
step. 

5. Main result 

5.1. Convergence Theorem 

We assume that the unique mild solution X of ADR problem ^ is the 
classical solution of ^ i.e. X is twice continuously different iable with respect 
to x and differentiable with respect to t. This assumption is necessary to 
achieve optimal orders of convergence in time and space. 

Theorem 5.1. Consider the mild solution X of the ADR M) and the nu- 



merical solution (53) given by the ETD1 scheme. Assume that Assump- 



tion ^2 is satisfied and the reaction function R satisfies Assumption 2A 



Set X h = PhX , assume that X Q e D(A), then the following estimate holds: 

\\X(t m )-X™\\ 0>h <C(At + h), 
where C = C(fl, X, R, D, q, T, Ci, (2)- 



The proof follows in Section 5.3, but we need some preparatory results 
first. 
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5.2. Preparatory results 

Proposition 5.2. [Interpolation error] 

Let T be an admissible mesh in the sense of Definition 3.1 and J~h its 
dual Delaunay triangulation ({x{\ are vertices of Th) ■ Let Ih '■ C(Q) — » Vh 
defined by 

I h {u) = ^\i( Xi V Xi , u e C(Tt) (54) 

ieT 

where {y2 x .} ig7 - is the nodal basis corresponding to {xj} ig7 - in the sense of 
finite element method ((p Xi (xj) = 8^). If u G C 2 (Q), then there exists a 
positive constant Cq(u) such that the following estimate holds 

\\u-I h (u)\\<C (u)h 2 . (55) 

Ifue C([0,T},C 2 (U)), then 

\\u(t) - I h (u(t))\\ <C (u,T)h 2 , VtG [0,T]. (56) 

Proof. See (27J Section 3.4, Theorem 3.29, page 139 or Exercise 3.25 page 
147] or [HJ Theorem 17.1, page 132]. □ 

Proposition 5.3. [Smoothing properties of the semi group[23]J 

Let ft > and < 7 < \, then there exists C > such that 

\\A^S{t)\\ L(L2m < Cr 13 for t>0, 
\\A-y{I-S(t))\\ L{L 2 m < CV for t>0. 

In addition, the following results hold 

A p 3(t) = S(t)A? on V(A^). 

If p > 7 then V(A 13 ) C V(A 1 ). 
\\D\S(t)v\\ p < cr 1 -^-^' 2 \\v\\ a , t > 0, v e V(A a ' 2 ) I = 0, 1, 

where D\:=^, \\.\\ a := \\A«' 2 .\\. 



Lemma 5.4. Let X be the mild solution of given in (13). Let t\,ti G 
[0, T], ti <t2- The following estimates hold: 
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(i) IfX G V(A) then 

\\X{t 2 ) - X{h)\\ < Cih-hY-U \\X \\ 2 + sup \\R(X(s) 



0<s<T 

for e G (0, 1/2) small enough. 



• (ii) If Xq G T>(A) and R satisfies the Lipschitz condition in (12) then 
\\X{t 2 )-X(ti)\\ < C(t 2 - h) (\\X \\ 2 + sup \\R{X(s)] 

\ 0<s<T 

Proof. Part (i). 

Consider the difference 

X(t 2 )-X(ti) 

= (S(t 2 ) - S'(ti)) X + ( [ 2 S(t 2 - s)R{X{s))ds - [ 1 5(ti - s)R(X(s))ds 



= I + 11, (57) 

so that \\X(t 2 ) - X^W < \\I\\ + ||77||. We estimate each of the terms 
and ||//||. For ||/||, using Proposition 5.3 yields 

||/|| = ||5 , (ti)A- 1 (I-5(* 2 -t 1 ))A 1 Xo|| < tf(t 2 -fi)||X || 2 . 
For the term II, we have 

// = [\s(t 2 -s)-S{t 1 -s))R(X(s))ds+ [ 2 S(t 2 - s)R(X(s))ds 

JO Jt! 
= Ih + Ih, 

with 

\\II\\<\\Ih\\ + \\II 2 \\. 
We now estimate each term \\Hi\\ and H-Z^H- For 

117/xll = || [\s(t 2 -s)-S(t 1 -s))R(X(s))ds\\ 
Jo 

< f 1 \\(S(t 2 -s)- S(tx - s))R(X(s))\\ds 



< I I' \\(S(t 2 -s)- S(h - s))\\ L(L , m ds) ( sup \\R(X(s) 

J \0<s<T 
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For e G (0, 1/2) small enough, using Proposition 5.3 yields 

< ( f 1 \\S(t 1 -s)A 1 ^A- 1+ ^l~S(t 2 ~t 1 ))\\ L{LHn)) ds) ( sup \\R(X(s) 

\Jo J \0<s<T 

< ( f 1 \\A 1 - e S(t 1 -s)A- 1+e (I-S(t 2 -t 1 ))\\ HL 2 m ds) ( sup \\R(X(s) 

\Jo J \0<s<T 

< C(t 2 - h) 1 -' ( f\ti- s)- 1+e ds] ( sup \\R(X(s) 

\Jo J \0<s<T 

< Cih-h) 1 -' ( sup \\R(X{s)) 

\0<s<T 

For || II 2 1|, using the fact that the semigroup is bounded, we have 
||//|| = || / 2 S(t 2 - s)R(X(s))ds\\ 

Jh 



< y* 2 \\S(t 2 -s)R(X(s))\\ds 

< (jf* \\R(X(s))\\ds 

< C{t 2 -h) ( sup \\R(X(s) 

\0<s<T 



Hence 

||//|| < HJ/xll + ||JJ 2 ||) ^Cfa-trf-'f sup (\\R{X(s))\ 

\0<s<T 

Combining previous estimations of ||J|| and \\II\\ ends the proof of part (i). 



Proof of part (ii) We consider again the difference in (57). The differ- 
ence with the proof of part (i) comes from the estimation of This time 
we write 

III = j\s{t 2 -s)-S{t 1 -s))R{X(s))ds 



o 



ti 

(S(t 2 - s) - Sfa - s)) (R(X(s)) - flpf(ti))) ds 



+ f\s(t 2 -s)-S{t 1 -s))R(X(t 1 ))ds 
Jo 



Hn + Ih 2 . 
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If R satisfies the Lipschitz condition given in (12), then using the result in 
part (i) together with Proposition 5.3 yields 



\ihi\\ < (J \\s(t 2 -s)-s(t 1 -s)\\ HL2m \\R(x(s))-R(x(t 1 )\\ds 

< C ( f 1 \\S(t 2 -s)- S(h - s)\\ L{L , m \\X(s) - Xih^di 



< C^tz-h) J (ti-s)- e tZa 

< C{t 2 -h), 

for e G (0, 1/2) small enough. We also have 

II//12II < ||iW*i)|||| [ t \s(t 2 -s)-S(t 1 -s))ds\\ L{L , m 

Jo 

< C\\ S{t 2 - s) - S(ti - s)ds || L{p{Q)y 
Jo 

Using the two transformations y = t 2 — s, y = t\ — s we find 
||//«|| = Cll/ 2 S(y)dy- T S(y)dy\\ L{L , m 

Jt 2 -t! JO 

fh rti i-t\ 

= Cll / S(y)dy+ / S(y)dy - j S(y)dy\\ L(L 2 m 
Jti-t\ Jti Jo 

rti rt%—t\ 

= c \\ I S(y)dy- / S(y)dy\\ L(L 2 m 
Jt x Jo 

< C{t 2 -h). 

The estimate of II\ ends the proof. □ 
Lemma 5.5. [Discrete Gronwall lemma fOfflj ] 

Let the sequence t n = nAt ^ T. If the sequence of nonnegative numbers e n 
satisfies the inequality 

71-1 

e n ^aAt^y-^ej + bt- ff , (58) 

for ^ (3, a < 1 and a,b ^ 0, then the following estimate holds: 

e n ^Cb r a (59) 
where the constant C depends on (3, a, a, T. 
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5. 3. Proof of Theorem 5. 1 



In the followng proof, and C[, i = 1, • • • , 6 are positive constants. 

Proof. We use the equivalence of the norms |.| and ||.||o,ft in Vh as the mesh 
T is regular [27] . 

Using the triangle inequality yields 

||X(i( m ) — X™||o,h < ||X(t m ) — Xh(t m )\\o,h + H^fe(^m) — -XjT||o,/i 

= I + 11. (60) 

Let us estimate /. Integrating the shifted version (by adding cqX in both 
size) of equation ^ over each control volume i E T and using the divergence 
theorem yields 

/x t (x,t)dx-^ / (DVX - qX) ■ n a da + c I Xdx. = R(x, t, X(x, £))dx.(61) 
For t G [0, T] , i G T and cr G using the same notation as in [7], let us set 



ri,a(t) -- 



mes^cr, 
1 

mes(cx) 



mes(cr)/i (T 



d(i,j 

[qi, a X i>+ - J a qX(t) ■ n a ] 
1 



(62) 



Pi (t) = X(xi,t) t- J.X(x,t)dx 



mesu 



mesu 



I i?(x, t, X(x, t))dx - R(xi,t, Xi(t)). 



Assuming that the unique solution X of ^ is the regular, Taylor expansion 
yields 

X t {x,t)=X t {x i ,t) + s i {x,t), \si{x,t)\ <C 1 (X,T)h 

f i X t (x,t)dx. = m.es(i)X t (x i ,t) + Si, Si = ^ Sj(x,t)dx, \Si\ < mes(i)Ci (X, T) h. 

Using Assumption 4.2[ mainly the fact that the coefficients of the diffusion 
tensor D and the Darcy's velocity q are regular, we can therefore used the 
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following similar results in [7] for general elliptic operators 
( \Ri^(t)\<C 2 (V,X,T)h, 

\r i}a (t)\<C!,(q,X,T)h, 
\RiAt)\ + M*)l < C 3 (q,D,X,T) h, 



(64) 



{ \ Pi (t)\<C' 3 (X,T)h. 

Using the fact that R is differentiable with respect to X and x, Taylor ex- 
pansion respectively yields 

mes(i)gi(t) 

i?(x, t, X(x, t))dx. - mes(z) J R(x i , t, Xi(t)) 



f(R(x, t, X(x, t)) - R{xi, t, Xi{t))) dx, 

</ i 

mes(i) (i?(x,, £, X(x i; £)) - Rfa, t, X^t))) + Jz^ t)(X(x, t) - X(x,, t))dx 



+ / Z 2 (x,i)(x - Xj)dx 

J i 

mes(z) (i2(xi, t, X( Xi , t)) - i2(xi, t, X^t))) + «(xi, t, X, i?), 



where 



^i(x,t) 



1 

3X 
- 1 <9i? 



x, + r(x - Xi), t, X(x i; t) + r(X(x, t) - X(x i; £))) 



Z 2 (x,*)= / — (x. i + r(x-x l ),t,X(x l ,t) + r(X(x^)-X(x i ,t)))dr 
K(x t ,t,X,R) = /Zi(x,t)(X(x,t) -X(x i ,t))dx + /z 2 (x,t)(x-x i )dx. 



If the solution X is differentiable with respect x, one more Taylor expansion 
yields 

Kxi,t,X,i?)| < mes(i)C 4 (i?,T,X)/i. 



Using the Lipschitz condition (12) yields 



mes(i)Qi(t) < mes(z) (C£(fi, i?, T, X) |X( Xi , t) - X(t) | + C 4 (i2, T, X) h) . (65) 
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Subtracting the first equation of (21) from (61) and using previous expres- 
sions yields 

mes (,)^l + £ G ia (t) + W^(t) + c mes(i)ei(t) 
at 



= J t (i*(x, t , X (x, t ) ) - R{xi , t, Xi (t) ) ) dx (66 ) 
+c mes(z)pi(t) - Yl mes(cr)(i2 i)(7 + r i)(T ) - S'j(t), Vz G T 



with 



f e i (t) = X(x i ,t)-X i (t), (€[0,31, 
G i>a (t) = -T a (ej(t) -e.it)), a = 
Gi, a (t) = r a e t (t), aeSiHdn, 
[ W il<T (t) = q iia (X(x a>+ ,t) - X a!+ (t)), 



(67) 



and 



Xj if q ■ n CT > 0, 
Xj if q ■ n CT < 0, 
Xj if qn a > 0, 



(68) 



a e £i n an. 



x CT , x CT G <9f2 if q ■ n CT < 0, 



Multipling equation (66) by e$(t) and summing for i G T yields 

mes(i) d(ej(t)) 



E 



+ E ei(t)(G i>a (t) + Wi, a (t)) + c mes(z)e?(£) 



Eei(t) (i?(x, t, X(x, *)) - i?( Xi , *, X,(t))) dx] 
ier 



(69) 



+E 



c mes(2)p i (t)e i (t) - ^ mes(cr)e i (t)(R i>a (t) +r it<T (t)) - e^S^t) 
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Let ef(t) a piecewise constant function defined by 

e r (t)=ei(t), for x6i, i G T, te[0,T\. 



As in the proof of Theorem 4.3, using the fact that 



C 5 <fM a <C 5 , We£, 



(70) 
(71) 



where \x a is defined in (17) and (16) yields 



f C 5 \\e T miT<\Mmi h <C' 5 \\er(t)\\lr 



\e T (t)\\l T =i:\D«e T (t)\* mes{(T) 



(72) 



where 



\D a e T (t)\ = \ei(t) - ej(t)\, if a = i\j, 

\D a e T {t)\ = \ei(t)\, if a G £i n dQ. 
Setting e CT]+ (t) = X(x CTj+ , t) — X Ut+ (t), as in the proof of Theorem 4.3 



(73) 



im- 



portant results from [7j yields 



Using (74) and (65) in the expression (69) yields 
( kmesW^ + \\e T (t)\\l h + c \\e T (t)\\l h < C>\\e T (t)\\l h + C 4 /* £mes(i)| ei (t)| 

+c C' 3 hJ2mes(€)\e i (t)\ + £ £ mes^e^i^) + r i)<r (*)) + d /i £mes(i)| ei (t)|. 

The continuity of the diffusion and advection flux at each interface yields 
-Ri.aOO = -Rj,a(t), r i)( r(t) = -r j)(r (t), for cr = z|j G £; ni . 
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Set 

R*(t) = \Ri,a(t)l r a (t) = \Ri,a(t)l «eT, cre£ int . 

Using the relation (64), the Cauchy-Schwarz inequality as in [7j for stationary 
elliptic problems and reordering the summation over the edges yields 

^^mes(a)e i (t)( J R i , CT (t) + r i)(T (t)) 
*C ^D ff e r (*)(i^(0+r ff (*)) 



Using the fact that ^mes(cr)c? (T ^ <i mes(f2)) and relation (72) yields 



^^mes^e^t)^^) + r ijCT (t)) 
< C 3 /i(nies(n)d)3||er(0lli ) r 

For any constant C > 0, Young's inequality yields 



(76) 



|C/i£mes(i)e^)| = | £(C ft mes(z)^)(mes(i)l ei (0)| < -||e r (t)||^ + -C 2 /i 2 mes(fi) 



ieT 



Ch\\e T 



ieT 



ih^ -C 2 /i 2 + -||e r 



2 

1,/r 



Using expression (j76j and (|77|) in expression (JT5J) yields 

<(^/2 + l)||er(t)||g, fc + C 6 /> 2 



EmesW ^Ml + || er(t) ||2 +Co || er 



2 

O./f 



(7S 



6 — ( - / 6l ( - / l) ^3) < - / 4; ^5 



Bounding the left hand side of expression (78) below yields 

.<*(<=?(»)) 



mesu 



ds 



^ (Ci + 2)||e r ( S )|| 2 + 2C 6 h 2 , \/s G [0,T]. (79) 
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Integrating both side of expression (79) through interval [0,t], < t < T 
yields 

\\e T (t)\\lh< \\er(0)\\l h + ^C 6 Th" + (C' A + 2) f \\e r (s)\\l h ds, Vt G [0, T].(80) 

Jo 



Applying the discrete Gronwall Lemma |5.5| yields 

IM*)||L < C(\\e T (0)\\l h + h 2 ) 



(81) 



C = C(fi,X,i?,D,q,T,Ci,C 2 ). 



Then 



/=||X(t m )-X fc (* m )||o, h = ||er(Ollo,fc<C'(||Xo-Xo fc ||o I h + /i). (82) 
If X h = PhX , we have 

\\Xq — Xoh\\a,h — \\Xo — PhXo\\o,h < Ch, (83) 
(see |27J), we therefore have 

/ = \\X(t m )-X h (t m )\\ 0A = ||e r (f m )llo,h < Ch. (84) 



Let us estimate II. From (50) and (52) we have 



771—1 



and 



k=0 



m—1 



tk+1 



X h (t m ) = S h (t m )X oh + J2 S h( f ™ - s ) P hR(^ X h {s))ds } (85) 



tk 



X™ = S h (t m )X oh + J2 S h(tm ~ 8)P h R(t k , X k h )ds 

i, n J th 



(86) 



k=0 Jtk 



The smoothing properties of the semigroup Sh in Proposition 5.3 and the 
equivalence |.| = ||.||o,h in Vh yields 

\\Xh(t m ) — X™\\ 0th = \\Xh(t m ) — X™\\ 



m—l 



< CJ2 I \\P h {R(s,X h (s))-R(t k ,X k h ))\\ds, 
k=o 
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and 



\\Xh(t m ) — X™\\o,h < 



m—l 

k=0 Jtk 

m-1 , , 

+£ 

fc=0 

//i + // 2 . 



(UPh^s,^))-^,^)))!!^ 



\P h (R(s,X(s))-R(t k ,X k h )) \\ds 



(87) 



For s G [0,T], using Proposition |5.2 yields 

||X,( S )-X( S )|| < ||X h ( S )-7 h (X( a )) + / fc (X( S ))-X( a )|| 

< (\\X h (s)-I h (X(s))\\ + \\I h (X(s))-X(s)\\) 

< (\\X h (s)-I h (X(s))\\+C(X,T)h 2 ). (88) 

Since Xh(s) — 7/j(X(s)) G 14, the equivalence |.|| = ||.||o,h and the uniform 
estimate of the term / in [0, T] yields 

\\X h (s)-I h (X(s))\\ = \\X h (s) - I h (X(s))\\ , h 

= \\X h (s) - X(s)\\ 0> h (by definition of \\.\\o,h) 
< C(ft,X,P,D,q,Ci,C 2 )/> (89) 

(90) 



Using the Lipschitz condition (12) with (89) and (88) yields 



Ih < C(n,X,R,B,q,T,(i,(2)h. 



We also have 

m— 1 

ih < E 



k=0 
m—l 



~'k + l 



+ £ 

k=0 Jtk 
= III + III 



\P h {R(s,X(s)) - R(t k ,X(t k )))\\ds 



P h (R{t k ,X(t k ))-R(t k ,X%)) \\ds 



(91) 



Using Lemma 5.4 and the Lipschitz condition (12) yields 



m—l 



tk+l 



\\X(s)-X(t k )\\ds 



k=0 Jtk 
m-1 „ tk+1 



< CJ2 / (s-t k )ds 



< C(T)At. 



(92) 
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Using Lipschitz condition (12), Proposition 5.2 and the equivalence ||.| 
IUIo,/i in V h 



m—1 



Hi < / " 7 * (*(**)) + h(X(t k )) - X k h \\ds (93) 

k=o 

m—1 ptk+i 

< CV / 1 \\X(t k )-I h (X(t k ))\\ + \\I h (X(t k ))-X k \\ , h ds(94) 

(95) 



fc=0 



< C(X,T)[h 2 + J2 \\X(t k ) - X*\\ Q)h ds 

k=0 ^ tk 



Then 



(m—1 ftk+i 
(At + h) + J2 1 \\X(t k )-X k \\ 0>h ds 
k=o , 



Combining estimates J and JJ yields 



\X(t r 



X™)\\o,h 



( m-l r t k+1 \ 

< C[At + h + J2 \Wt k ) - X k \\ 0>h ds . 

V k=0 Jt * J 



(96) 



Applying the discrete Gronwall Lemma 5.5 in (96) ends the proof 



□ 



6. Simulations 

The Key point in (ETD1) scheme is the computing of the exponential 
function, the so called ifi. It is well known that a standard Pade approx- 
imation for a matrix exponential is not an efficient method for large scale 
problems [29]. Here we use the real fast Leja points and the Krylov sub- 
space techniques to evaluate the action of the exponential matrix function 
(pi(—At Ah) on a vector v, instead of computing the full exponential function 
<fi(— At Ah) as in a standard Pade approximation. The details of the Krylov 
subspace technique are given in [181 El El HE] while more information about 
the real fast Leja points can be found in [201 EH El El 119] - 
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6.1. Space convergence 

The aims here is to analyse the convergence in space by comparing the 
numerical solution to the exact solution for a linear problem. We consider 
here a linear reaction diffusion with exact solution given in [15]. The domain 
is defined as Q = [L ,Li) x [L ,Li), L = 0.01, L\ = 2. The initial time 
is given as to = 0.01. This is necessary because the exact solution is not 
defined at the origin and at t — 0. The dispersion tensor D is diagonal and 
heterogeneous with coefficients given by 

D lt i(x,y) = D ulx 2 (x,y) efi 
D 2 ,2(x, y) = DquIv 2 (x,y) efi. 



A curved the velocity field (Figure 2(a)) is given explicitly by 



q = (Qx, q y ) T _ 

q x (x,y)=u x (x,y) eQ_ (97) 
Qy(x,y) = -u y (x,y) Gfi 

where Do = 0.1 and uq = 2. The finite volume mesh T is the rectangular grid. 
Initial and boundary conditions are taken according to the exact solution [15] , 
assuming an instantaneous release at a point (xo,yo), xo = 1.5, yo = 1.5. We 
take a fixed time-step of At = 1/3000. We use here the real the fast Leja 
to compute the action of the exponential function (pi. The final time here is 
T = 1. 



Figure 1(b) shows a convergence of order O(h) for the spatial discretisa- 



tion, which agrees with the predicted order in Theorem 5.1 



6.2. Time Convergence 

We use the upper 40 layers of the highly heterogeneous SPE10 case, which 
represents a braided fluvial North Sea oil field [1] for porosity and perme- 
ability fields. 

The domain is Q = [0, Li] x [0, L 2 ] x [0, L 3 ], the finite volume mesh T is 
the set of parallelepipeds with space discretization stepsize Ax = 20 ft, Ay = 
10 ft and Az = 2 ft. Dimensions are L x = 1200/t, L 2 = 2200 ft, L 3 = 40 ft. 
We use for the absolute tolerance error tol = 10~ 6 and m = 8 for the Krylov 
subspace dimension. 
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^ log( Space stepsize) 

(a) (b) 

Figure 1: (a) shows the streamlines of the velocity field, (b) shows the L 2 norm of the 
error at T — 1 as a function of the grid size h for linear advection diffusion. 

To have the velocity field we solve with our finite volume scheme the 
system 

• Vq=0, 

where the dynamical viscosity is fi = 0.3cp and the diagonal permeability 
tensor K comes from [TJ. We use the Dirichlet boundary at 

Y D = {{0} x {0} x [0,L 3 ]} U {{LJ x {L 2 } x [0,L 3 ]} , 

and homogenous Neumann boundary elsewhere such that 

( 3998.96 psi in {0} x {0} x [0, L 3 ] 
V = I 

{ 7997.92 psi in {L x } x {L 2 } x [0,L 3 ] 

-KVp(x,t) ■ n = in T N = dQ\T D . 

For the nonlinear advection diffusion reaction we use the following boundary 



28 



conditions 



X = in {{0} x {0} x [0,L 3 ]} x [0,7] 
X = 1 in {{LJ x {L 2 } x [0,L 3 ]} x [0,T] 

-(DVI)(x,t) -n = in rjvx[0,T] 

Xo = in Q (initial solution) 

where n is the unit outward normal vector to T^- The velocity field ob- 



tained in (98) is piecewise constant in the mesh T. For reaction function, 

— Xf3X 

we use the classical Langmuir sorption isotherm given by R{X) = — — , 

1 + XX 

with A = 1, = 10~ 3 . This reaction function is obviously Lipschitz. The 
reference or 'exact' solution is the numerical solution with the smaller time 
step size (At = 0.25). The final time here is T = 4096 seconds. We take a 



uniform dispersion(diffusion) tensor of D = 10 6 x I 3 . Figure 2(b) shows a 



convergence of order O(At) for the time discretisation, which agrees with the 



predicted order in Theorem 5.1 even we don't have any information about 
the regularity of the exact solution. In fact to obtain the optimal order of 
convergence in time, the exact solution should be differentiable respect to 
the time, which may be the case here since the initial solution is infinitely 
differentiable and we are dealing with parabolic problem. 
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